18/05/24
Version with polarization dependence.
22/03/24
Alignment frame calcs preliminary:
Note all plots are interactive, mouse-over for data points, and use toolbars for zoom/selection, and widgets for browsing data types (where applicable).
For general methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
Now with ADM class.
from epsproc.classes.alignment import ADM
from pathlib import Path
dataPath = Path('~/ePS/N2O/alignment/Wigner D polyatomic normalized')
ADMtype = '.txt'
ADMin = ADM(fileBase=dataPath.expanduser(), ext = ADMtype)
ADMin.loadData(keyType='setADMsT', normType='wignerDpoly')
# import epsproc as ep
# ep.setADMs
ADMin.plot(width=1000)
# ADMin.plot(xlim=(30,50), width=800) # Auto-stack plus pass args & zoom in on specific region (note slider will reset region with interactive zoom)
# Plotters
# from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
import warnings
# warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
# # Scan for subdirs, based on existing routine in getFiles()
# fileBase = Path('/home/paul/ePS/OCS/OCS_survey') # Data dir on Stimpy
fileBase = Path('~/fock-mount/globalhome/eps/N2O/N2O_valence').expanduser() # Data dir on Jake
# TODO: fix orb label here, currently relies on (different) fixed format
data = ePSmultiJob(fileBase, verbose = 0)
data.scanFiles()
# data.jobsSummary()
Using EfiledPol class, set various ellipticities and fields $(E_x,E_y)$. For MF calcs, compute for all polarization geometries corresponding to fields aligned along $(x,y,z)$ molecular axes.
Note default case assumes $z$-axis alignment, and $E_y$ for the field plots here is reoriented along the alignment axis in the AF calcs. See E-field pol docs for details and control over geometries.
18/05/24: debugging for alignment... in this case RX is not used in calcs, so need to ROTATE FIELD TERMS prior to PAD calcs.
from epsproc.efield.epol import EfieldPol
import numpy as np
import matplotlib.pyplot as plt
# Mulitple ellipticities and rotations
# For fields, set with format = [amplitude,azimuth (0,pi), ellipticity (-pi/4,pi/4)]
states = 3
# angles = np.linspace(0, np.pi, states) # For range of angles
angles = np.zeros(states) # no rotation
maxEll = 0.5 # Set max ellipticity to use (0 - 1)
ellipticities = np.linspace(0,maxEll*np.pi/4, states) # Set for 0 - max ellipticity
ell = np.c_[np.ones(states),angles,ellipticities]
labels = (ellipticities/(np.pi/4) * 100).round(2) # Set labels at %age ellipticity
# Set fields
Eell_multi = EfieldPol(ell=ell, labels = labels)
# Plot fields
Eell_multi.plot(figsize=(6,6), draw_arrow=False)
plt.legend(labels, loc='upper right')
# Optional - set field orientation for PAD calcs. If not set defaults to Ey>z-axis rotation.
# Eell_multi.setOrientation
Eell_multi.setOrientation()
import epsproc as ep
ep.setPolGeoms
from epsproc.geomFunc.geomUtils import EfieldPolConfig
ep, p, RX, EPRX, EfieldPol = EfieldPolConfig(Eell_multi)
ep
p
RX
EPRX
Eell_multi.YLM
%matplotlib inline
# Note some build envs may need mpl inline reset here!
# Plot light fields - input
# Eell_multi.plotSph(dataType='rot', facetDim=['t','Euler'])
from epsproc.sphPlot import sphFromBLMPlot
Itp, fig = sphFromBLMPlot(Eell_multi.YLMrot.sel(t=0).squeeze(drop=True), plotFlag = True,
backend='mpl', facetDim='Euler')
Itp, fig = sphFromBLMPlot(Eell_multi.YLMrot.sel(t=0).squeeze(drop=True), plotFlag = True,
backend='mpl', facetDim='Euler')
Eell_multi.YLMrot
Eell_multi.EPRX
from epsproc.sphCalc import setPolGeoms, setBLMs, TKQarrayRotX
TKQarrayRotX(EPRX, RX)
EPRX
break
TO DO May 2024: E-FIELD ROTATION... calcs below seem to be for xpol case otherwise (XS show opposite alignment dependence to previous calcs)
Eell_multi.YLM
Eell_multi.YLMrot
Eell_multi.YLMrot.swap_dims({"Euler":"Labels"}).unstack('BLM').sel(Labels='z')
# Test EPRX setting from rotated case - currently missing...
from copy import deepcopy
Etest = deepcopy(Eell_multi)
# Test replacing epDict with rotation (subselection)
subselection = Eell_multi.YLMrot.swap_dims({"Euler":"Labels"}).unstack('BLM').sel(Labels='z').squeeze()
Etest.epDict = {'ep':subselection.values,
'p':subselection.m.values,
'RX':Etest.RX}
Etest.epDict
Etest.calcEPR()
Etest.epDict
hasattr(Etest,'YLMrot')
Etest.EPRX.unstack()
Eell_multi.epDict
# Set ADMs to use for calc - set for specific T and downsample
# Subselection, note temperature is set by dataKey here
temp = '10K'
trange = [38,42]
# ADMin.subsetADMs(dataKey = temp, trange=[35,45],tStep = 5)
ADMin.subsetADMs(dataKey = temp, trange=trange,tStep = 5)
# Plot subselection
ADMin.plot(keys='ADM')
# Set to main data structure for calcs
data.data['ADM'] = ADMin.data['ADM']
tpoints= data.data['ADM']['ADM'].t.size
# Basic routine (no grouping)
thres = 1e-5
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-2) # OK for 41 t-points (x2 downsample), but missing some orbs and t-points in output!
# data.AFBLM(AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=thres) # OK for 41 t-points (x2 downsample), still some missing regions for thres=1e-2
# mostly contiguous for thres = 1e-3,1e-4
# EfieldInput = Eell_multi
EfieldInput = Etest # Testing frame rotation...
# With pol...
data.AFBLM(AKQS = data.data['ADM']['ADM'].real, EfieldPol = EfieldInput,
selDims = {'Type':'L'}, thres=thres)
# For thres=1e-3, 3 pol states, 41 t-points, requires ~ 50Gb RAM/10 mins. Some dicontinuities
# May need grouping for larger calcs.
# Thres = 1e-4, 3 pol states, 16 t-points, similar requirements
# Testing 22/05/24 - implementation of field frame rotation for EPRX
# Save full class with Pickle
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
import pickle
outfile = f'N2O_AFBLM_polDep_{temp}_t{trange[0]}-{trange[1]}_tp{tpoints}_{thres}_frameRot_220524.pickle'
with open(outfile, 'wb') as handle:
# data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = [] # UGH, this breaks pickle unless removed.
pickle.dump(data.data, handle, protocol=pickle.HIGHEST_PROTOCOL)
print(f'Saved data.data with pickle, file: {outfile}')
# # Load data
# thres = 1e-5
# # infile = f'N2O_AFBLM_{temp}_{thres}_220323.pickle'
# # infile = 'N2O_AFBLM_10K_0.0001_220323.pickle'
# infile = 'N2O_AFBLM_10K_1e-05_220323.pickle'
# import pickle
# with open(infile, 'rb') as handle:
# dataIn = pickle.load(handle)
# data.data = dataIn
NOTE: interactive plots. Roll over plots for values and tool options. For heatmaps use the "box select" tool on the histogram at the side to change the colour mapping scale.
Remarks:
# v2 with pkg code...
# hvTest = data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = False, addADMs = False, XS=True, renderPlot=False, returnPlot=True)
data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = False, XS=True, clim=(-1,2))
# Plot some results
# Example with post-plotter subselection for faster plotting...
# Create HV object, but don't plot
data.BLMplot(dataType='AFBLM', xDim='pol', thres = 1e-2, backend='hv', width=700, ylim=(-1.5, 2.5),
renderPlot=False)
# Replot from Holoviews dataset with selectors
# Note this may need hv or hvPlotters to be loaded
from epsproc.plot import hvPlotters
xDim = 't'
# data.plots['BLMplot']['hvDS'].select(Eke=[1.1,10.1,20.1]).select(Orb='orb9_P').to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm','pol'])
# data.plots['BLMplot']['hvDS'].select(Eke=[1.1,10.1,20.1]).to(hvPlotters.hv.Curve, kdims=xDim).overlay(['l', 'm'])
data.plots['BLMplot']['hvDS'].to(hvPlotters.hv.Curve, kdims=xDim).opts(ylim=(-1.5,2.5), width=1200).overlay(['l', 'm'])
data.data['ADM']['ADM']
# from epsproc.plot import hvPlotters
# hvPlotters.setPlotDefaults(fSize = [700,300], imgSize = 600)
data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = True, XS=False,
clim=(-1,2), cmap='vlag')
data.plots['BLMplot']['XR'].attrs['harmonics']['dimList']
def testKWARGS(**kwargs):
if 'clim' in kwargs.keys():
print(f"clim = {kwargs['clim']}")
else:
print("Nope")
testKWARGS(nothing='thing')
testKWARGS(clim='thing')
testKWARGS(clim=(-1,2))
plotObj = data.BLMplot(backend='hv', xDim = 't', hvType='heatmap', addHist = True, addADMs = True, XS=False,
clim=(-1,2), cmap='vlag', returnPlot = True, renderPlot = False)
plotObj.hist(bin_range=(-2, 2))
# bin_range=(-5, 5)
break
# hvTest.hist()
# hvTest.overlay(data.plots['BLMplot']['XR'].attrs['harmonics']['dimList']) #.opts(**kwargs)
abs
subset = data.plots['BLMplot']['XR'].copy()
# selDims =
import epsproc as ep
ep.matEleSelector(subset.unstack(), thres=thres, inds = {'l':slice(2,20)}, drop=True) #, dims = contiguousDims, sq = sqSelector)
# ep.clea
subset.where(subset.l>0, drop=True)
from epsproc.plot import hvPlotters
import numpy as np
# Set data - use existing plot routine but skip line outs (quite slow, replotted at end of notebook)
data.BLMplot(backend='hv', xDim = 't', renderPlot=False)
# Heatmaps for cross-sections
xrAF = data.plots['BLMplot']['XR'].copy()
Etype = 'Eke'
# Set opts to match sizes - should be able to link plots or set gridspace to handle this?
hvPlotters.setPlotDefaults(fSize = [750,300], imgSize = 600)
cmap = 'coolwarm' # See http://holoviews.org/user_guide/Colormaps.html
# Get cross-sections
XS = np.abs(xrAF.XSrescaled) # Should be in MB
XS.name = 'XS/MB'
hvXS = hvPlotters.hv.Dataset(XS.fillna(0))
# Plot heatmap + ADMs
(hvXS.to(hvPlotters.hv.HeatMap, kdims=['t',Etype]).opts(cmap=cmap).hist() + data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')).cols(1)
# Heatmap for l=2,4,6
# NOTE - clipped cmap to physical values with clim=(-1,2), but some spikes.
# BLM = np.abs(xrAF.XSrescaled) # Should be in MB
# XS.name = 'XS/MB'
# Filter out small XS regions?
xrAFfiltered = xrAF.where(np.abs(xrAF.XSrescaled) > 0.01)
hvDS = hvPlotters.hv.Dataset(xrAFfiltered.unstack('BLM').fillna(0).squeeze(drop=True))
# hvDS = hvPlotters.hv.Dataset(xrAF.unstack('BLM').fillna(0).squeeze(drop=True))
# Plot heatmap + ADMs
(hvDS.select(l=[2,4,6], m=0).to(hvPlotters.hv.HeatMap, kdims=['t',Etype]).opts(cmap=cmap, clim=(-0.5,1)).hist() + data.data['ADM']['ADM'].unstack().squeeze().real.hvplot.line(x='t').overlay('K')).cols(1) #.overlay(['l','t']).layout('Orb').cols(1)
# Line-outs
# Note this is quite slow to render
data.BLMplot(backend='hv', xDim = 't', ylim=(-1,2), filterXS=0.01)
hvDS.select(l=[2,4,6], m=0)
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews', 'pandas'])
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc
22/03/24
v1, quick run based on previous OCS code plus a few recent function updates.
TODO:
break
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))
#*************** Throw everything in method - breaks for large datasets (RAM issues)
# data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}) # NOTE THIS USES independently set ADMs!
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3) # OK for 61 t-points (x2 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-3, thresDims = ['Eke','t']) # NEED TO ADD DIM CHECKING TO thresDims
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-4) # 1e-4 working OK for 30 t-points (x4 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-5) # 1e-5 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=1e-6) # 1e-6 working OK for 16 t-points (x8 downsample)
# data.AFBLM(keys = keys, AKQS = data.data['ADM']['ADM'].real, selDims = {'Type':'L'}, thres=None) # FAILS for 16 t-points (x8 downsample) at 60Gb RAM on Jake (Swap file size issue????)
# OK for 8 t-points (x16 downsample) at ~55Gb RAM peak on Jake.
# thres=1e-4) BREAKS KERNEL #, thresDims=['Eke','t']) # TODO: currently thersDims must be in matE only for AFBLM routine.
# 21/07/22 version - above only gives L=2 max, missing L=4,6. But bumping threshold to 1e-4 kills calculation...
# data.AFBLM(keys = keys, selDims = {'Type':'L'}) # Default case... sets/uses ADMs per key
# UPDATE: testing with scale factor for ADMs, and forcing to real - might be better?
#*********** Using grouping to avoid memory issues
# QUITE slow with current code, but more robust
# TODO: make this better with basis set recycle (cf. fitting code)
#
import pickle
def calcAFBLMGroup(data,ADMs,keys,AFBLMdict):
"""Wrap for XR groupby"""
# Class routine - returns to main datastruture
data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'}, thres=None) #.map(standardize)
# Copy to new structure & restack
# keys = data._keysCheck(keys)
# AFBLMdict = {}
# Restack data to lists for merge
for key in keys:
AFBLMdict[key]['AFBLMlist'].append(data.data[key]['AFBLM'].copy())
return AFBLMdict
AFBLMdict = {k:{'AFBLMlist':[]} for k in keys}
# tBins = np.arange(38,39.1,0.5) # Test range
tBins = np.arange(trange[0], trange[1] + 0.1, 0.5)
# trange
import time
for item in data.data['ADM']['ADM'].groupby_bins("t", tBins):
start = time.time()
print(f'Calculating for group: {item[0]}, {item[1].t.size} t-points.')
AFBLMdict = calcAFBLMGroup(data, item[1], keys, AFBLMdict)
end = time.time()
print(f'Delta: {end - start}s')
# Concat XR lists
for key in keys:
AFBLMdict[key]['full'] = xr.concat(AFBLMdict[key]['AFBLMlist'],"t")
AFBLMdict[key]['full'].name = f'{key}_AFBLM'
# Push back to main class
data.data[key]['AFBLM'] = AFBLMdict[key]['full'].copy()
# Save object with Pickle
with open(f'OCS_{key}_AFBLM_220722.pickle', 'wb') as handle:
print(f'Saving {key} with pickle')
pickle.dump(AFBLMdict[key]['full'], handle, protocol=pickle.HIGHEST_PROTOCOL)
# Save full class with Pickle
# NOW GIVING ERRORS ON RERUN (although worked previously), AttributeError: Can't pickle local object 'plotTypeSelector.<locals>.<lambda>'
# NOT SURE WHERE THIS IS COMING FROM AS YET
# UPDATE: IN data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails']
# TODO: GENERAL DEBUG FOR ACCIDENTAL FUNCTION PASSING HERE!!!!
# import pickle
with open(f'OCS_AFBLM_220722.pickle', 'wb') as handle:
data.data['ADM']['plots']['ADM']['pData'].attrs['pTypeDetails'] = [] # UGH, this breaks pickle unless removed.
pickle.dump(data, handle, protocol=pickle.HIGHEST_PROTOCOL)
print(f'Saved data (full class) with pickle')